
SOUTHWEST RESEARCH INSTITUTE 


SAN ANTONIO, TEXAS 



SOUTHWEST RESEARCH INSTITUTE 
8500 Culebra Road, San Antonio 6, Texas 


LIQUID SLOSHING IN A SPHERICAL TANK FILLED 
TO AN ARBITRARY DEPTH 


by 


Wen-Hwa Chu 


Technical Report No. 4 
Contract No. NAS8-1555 
SwRI Project No. 6-1072-2 


Prepared for 

National Aeronautics and Space Administration 
George C. Marshall Space Flight Center 
Huntsville, Alabama 


15 December 1962 

APPROVED: 



H. Norman Abramson, Director 
Department of Mechanical Sciences 



TABLE OF CONTENTS 

Page 

FORWARD in 

SUMMARY 1 

INTRODUCTION 2 

I. MATHEMATICAL FORMULATION 6 

A. Kernel Function 6 

B. Eigen Function 8 

C. Sloshing Force and Pressure in Translational 

Oscillation 9 

D. The Moment Under Translational Oscillation 11 

E e Pitching Oscillation 11 

II. NUMERICAL METHOD 

A. Approximate Determination of the Kernel Function 

at a Finite Number of Points 14 

B. Determination of Eigen Vectors, f and Natural 

Frequencies 19 

C. Evaluation of Force 2 0 

D. Precision Problem 21 

III. EXAMPLE: FUEL SLOSHING IN A QUARTER -FULL 

TANK 24 

IV. CONCLUSIONS AND DISCUSSIONS 27 

ACKNOWLEDGEMENT 28 

REFERENCES 29 

NOMENC LATURE 3 1 

APPENDICES 34 

Appendix I. Analytic Expression for HoCf'f*) 34 

Appendix II. Analytic Expres sion for 43 

Appendix III. Subroutines DKEF and NEFF 

(WIZ Program) 50 

Appendix IV. Derivation of 7 7 n 84 

Appendix V. X-Force Acting on the Tank by 

Integration of Pressure 58 


TABLE OF CONTENTS (Continued) 

Page 

FIGURES 60 

1 Graphical Illustration of Some Nomenclatures 60 

2 Surface of Integration 61 

3a Moment About Axis of Rotation 62 

3b Motion of Tank in Pitching (Rocking) Oscillation 63 

4a Comparison of First Natural Frequency with Data 

by Abramson, et al. , (Ref. 15) 64 

4b Comparison of the First Three Natural Frequencies 

with Experiments by Stofan- Armstead (Ref. 10) 65 

5 Comparison of Force Response for Quarter-Full Tank 

with Experiments by Abramson, et al. , (Ref. 15) 66 

6 Natural Frequencies Given by Riley-Trembath (Ref. 17) 67 


ii 


FORWARD 


The work presented in this report is the culmination of efforts 
on the part of the author that have extended over a period of 
two years, with partial support being received from several 
sources. The analysis was originally undertaken as part of a 
theoretical study under Contract DA-23-072-ORD- 1251 spon- 
sored by the Army Ballistic Missile Agency. The work was 
continued, and largely completed, under the program of inter- 
nal research supported by Southwest Research Institute {Pro- 
ject 1059-2). Because of its relevance to experimental work 
already completed under the present program (NASA-MSFC 
Contract No. NAS8-1555) and published as Technical Report 
No. 2, it was felt highly desirable that this work should also 
be issued and distributed under auspices of this contract. 

The material in this paper is also presented in Part I, Section 3, 
of a dissertation titled "Some Contributions to Unsteady Hydro- 
dynamics in Engineering, " submitted to the graduate faculty of 
The Johns Hopkins University in partial fulfillment for the 
requirements for the degree of Doctor of Philosophy. 
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LIQUID SLOSHING IN A SPHERICAL TANK FILLED 
TO AN ARBITRARY DEPTH 

SUMMARY 


The kernel function for liquid sloshing in a spherical tank 
filled to an arbitrary depth is shown to be related to the Green’s 
function of the second kind and is constructed successfully by numer- 
ical means. Natural frequencies are then computed as eigen values 
of a matrix. Eigen functions are obtained at a finite number of points 
as the eigen vectors which are sufficient for approximate evaluation 
of the force acting on the container. Simple formulas of force and 
moment are given for both pitching and translational oscillation 
under a fixed gravitational field. Finally, comparisons of predicted 
natural frequencies and force response with experiments for a 
quarter -full tank are also given. 



INTRODUCTION 


Disturbances on a rocket or missile can induce sloshing of 
fuel in a partially filled tank. It in turn exerts excitation forces on 
the vehicle and in some cases can be detrimental to the trajectory or 
even results in loss of control. Sloshing in a circular cylindrical tank 
has been widely investigated with and without damping. To facilitate 
dynamic analysis, an equivalent mechanical model for circular tank 
is given in Reference 2. For a spherical tank, an ingenious semi- 
numerical method was given in Reference 1. However, the problem 
is only solved for three special cases, namely, nearly full, nearly 
empty, and half-full tanks. The restriction is due to the lack of the 
Green* s function of the second kind (Neumann function) for the spher- 
ical bowl. Although the Green* s function of the first kind for the 
spherical bowl is given in Reference 3, it is doubtful that a simple 
expression for the Green* s function of the second kind exists in the 
toroidal coordinates, since the normal derivative on the spherical 
cap is a combination of two derivatives in this coordinate system. 

The sequence method given in Reference 4 is convergent for Green’s 
function of the first kind but may diverge for the second kind. One 
may resort to Laouville -Neumann method (series method, Ref, 5) and 
prove it converges. But when the Green’s function on the boundary 
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is desired, the kernel function is singular; thus it becomes increas- 
ingly more difficult to evaluate when more terms are needed. If we 
do not employ the Neumann function, an integral equation on the free 
surface is also obtained. Unfortunately, the eigen functions no longer 
satisfy the necessary orthogonal relationship (Ref, 6), thus they are 
the desired eigen functions only if the Neumann function is employed 
(Ref, 1), In this paper, a numerical scheme is devised to determine 
the desired kernel function, which is one component of the Neumann 
function, and then apply the same procedure as given in Reference 1 
to evaluate the sloshing characteristics. Considerably more work is 
required to calculate the pressure on the wall, although in principle 
this can be done. 

After the theory in the present paper was developed, some 
other approaches have been published. One approach (Ref. 7) seeks 
the variational solution based on Hamilton’s principle through Rayleigh - 
Ritz method*. Since only an integrated free surface condition was im- 
posed, it is somewhat doubtful that accurate prediction of force response 
or pressure can be assured (Ref. 8), although error in the lowest mode 
frequency was less than one per cent for a flat cylindrical tank. In 
another approach (Ref. 9) finite difference techniques were employed to 

This method has been applied to spherical tank by Riley and Trembath 

whose results are shown in Figure 6. 
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seek eigen values in a boundary condition by three different methods . 
Method I and Method III (Ref. 9) use either Rayleigh quotient or Rayleigh - 
Ritz procedure, but are somewhat inferior (Ref. 9) to the Rayleigh-Ritz 
procedure applied to the continuous domain. Method II (Ref. 9) converts 
the problem into an equivalent matrix eigen value problem by eliminating 
the points outside the free surface through an inversion of matrix if the 
number of the other points is small, or through an influence coefficient 
type calculation if otherwise. In the latter case if there are N points 
on the free surface, N boundary value problems should be first solved 
(say by successive over -relaxation) before reduction to the eigen value 
problem of a N x N matrix. Depending to a large extent on the number 
of net points required for a desired accuracy (say, 3 figures in frequen- 
cies and force response), the computing time (based on estimation on 
a GE 225 computer)* of the last method for a spherical tank seems to 
be comparable to the present method. On the other hand, although 
further (significant) acceleration of the rate of convergence of the sub- 
routines in the present method in the present problem may be quite 
difficult, an alternative numerical scheme devised is expected to reduce 

* T • 

It is estimated under the assumption that there are 20 free surface 
points and 300 total net points with 120 iterations for each boundary 
value problem (based on experience of a similar problem) and 
average speed for 5 multiplications, 4 additions, and one additional 
multiplication or division at each point in each iteration. There are 
other estimates based on experiences which yield approximately the 
same magnitude of computing time. 
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the computing time to one -half or further. Finally, Reference 10 has 
also been published in which the kernel function is constructed empiri- 
cally, based on knowledge for half -full and full tank. 

The purpose of the present paper is mainly to predict the natural 
frequencies and force response and to show how kernel functions are re- 
lated to the Neumann function on the boundary and can be constructed 
numerically for a spherical tank. Analogous extension to other con- 
figurations or other problems may be possible but will not be treated 


in this paper. 


MATHEMATICAL FORMULATION 


Kernel Function 


Let G(P_Q) and (j e (P,Q ) be the Green's function of the second 
kind for the interior of the given spherical bowl (Fig. 1) and the sphere, 
respectively: (a) Both G(P, Q ) and Q) possess continuous second 

derivatives and satisfy the Laplace equation inside the bowl and the 
sphere, respectively, except the point P—Q ; (b) Both G and G o possess 
a unit sink, — at P* Q inside the bowl; (c) ~~ 4nO? ° n 




whole surface of the sphere, =* ~ on the surface of the 

bowl, R and F; (d) be that given in Reference 7; Gr(.Pj Q) satisfies 


the normalizing condition J Q) c/S r — O (Ref. 11). Following 

ft+F 

these conditions, it is well known (Ref. 11) that the Neumann function G 
is symmetric as well as G a , i. e. , G(P,Q) *=G(Gl,P) G a (P.Q)= GJQ,P) 
When Q are both interior points, analogous to the proof of symmetric 
properties, one has 

G(Q.p>- 6.<.P.Q) - { 9.(10) - 9(1. P) J J Sl 

- *, G.(ro)^-/ r sn, p; 

~b F riy 4 

which is an integral equation governing <5 ( P, Q) where P, Q is inside 
the bowl, not on F and R. 


For values of the Green's function with p J Q (P^Q) both on F, 
not on R, apply directly the divergence theorem to the surface shown 
in Figure 2. Since there is an infinitesimal semi-sphere around the 
sinks at P and Q respectively, one finds 

j G( Q,P) -£■ <7, (p, Q) = f [G.d.Q^-GfcP) ^J7^] dS x 

R+F 

P*i*Q 

By making P and Q in Equation Cl] approach P and Q on the free sur- 
face along its normal, Equation Cl] can be reduced to Equation [2J. 

In Reference 1, for fuel sloshing in a spherical tank, only those 
eigen functions proportional to cos 9 are needed: one shall see in the 
next section that it is sufficient to know one component H(P'Gt) of the 
Neumann function G(P, Q) to determine the sloshing characteristics. 
Let 

Z7t 27T 

H(p, O) = 4 f 0 J 9 G(P,Q) cos 9p c 03% </9 p de & 
zn 

K(P, * £ jf J o «« Gp cos J& p 


k'(p,Q)= 


nr 



G(PQ)-G e (P,Q)] cos&p COS 6^ dQpde ^ = H(p. Q)-H 0 (P, Q) 


Since G and G 0 are symmetric functions, H , H v and thus ^ are 


symmetric functions. 


For points P, Q corresponding to P and Q respectively, inside 


the spherical bowl, Equation [l]j can be integrated to yield 

/i,(p,Q) = -/_ {%(?.&) A, 

F F 

for which the reversing of orders of integration are applied and can be 
justified by carrying out the details. The function F is defined by 

zn r \ z1 r / T 

fJlQ) !— f ±f( cosc^i) . 

cos Sri « « *7 J Tn x 




which is a nonsymmetric function as 

Similarly, if both P and Q are on F, not on R, integration of 
Equation [2} yields 

ua? a z(/t.fD] [j/qi 

[ d/i fla /q)] 


where and // # are given in Appendices I and II. For = 0 , 

— 0 , almost everywhere on F, hence for half- sphere h,(R. P)= (P, //=/4 

which is in agreement with Reference 1, 


B. Eigen Functions 

The eigen functions <j> n are assumed to possess the following 
properties: (a) ^ is regular inside the bowl and 7^ =• o 

(b) T^= 0 onR ’ 1 7%T=X, ( k( I ) on F, (c) <fc(p) = %(P) cose p (Ref. 1). 

The last condition is appropriate for translational oscillation of the tank. 
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Analogous to Equations fl] and (Zl 

fn(P) t&dSt 

when P is inside the spherical bowl, R+ F. 

i &(p)= f g(i, p) Z id) </s 7 

when P is on F. 

Analogous to Equations [4] and £5^ by integration 

%(?>-( Hfr & && , p not on F. 

J F 

t> _ _ 

Jf' ’ p on F - 

This shows only H(P,Q) is needed for the pertinent eigen functions. 


C. Sloshing Force and Pressure in Translational Oscillation 

By introducing a displacement potential relative to the tank 
the sloshing force acting on the container is derived 
from the Lagrangian* s equation in Reference 1, namely 


= 2 fa d„ 

where * mass of the liquid = £ 3* £ £a *+ 30 % -Z*] 

* = fl m £ dFl 

“» ” J F ds = Vk 'J, t<f,) tl "ft 


[73 

C8] 

f 9 ] 
flO] 

tH ] 

Ola] 

(lib] 
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Q.M — 


s_._g» 

tO/> **ri 


u 


tO* 


-/ 


The velocity potential 


[lie] 


t * * Q»(*)h(r.'ke) +UX [12] 

The pressure on the container 

-f> 2 Q„(t)<k(a,<p,e)- f>u* ri3] 

within the accuracy of the linearized theory. Equation [11] can also 
be obtained directly by inte gration of pressure (Appendix V). 

Once ^ on F is evaluated, one may employ Q 0 to obtain 
^(P) from 




COS 0 p 


fjp) « f_ d* (?,P) X. -f_ % (1 J P) %(x)ds x 

F F 

d <f> 

The integral on R dropped out as — 0 on R and 
on R. 


2G 0 

— ~ = constant 


For P on R, not on F, the integrands of the integrals in 
Equation [l4bj are nonsingular, hence f n (p) can be calculated by 
well-known numerical methods. For contact points both on R and F, 
the value of may be obtained by evaluation of the integral by 


[l4a] 

[14b] 


midpoint formula. 
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D. The Moment Under Translational Oscillation 

For translational oscillations, the velocity potential is propor- 
tional to cos & and the flow is antisymmetric. It produces a horizon- 
tal force in the x-direction and a couple about the center of 

the tank (Fig. 3a). There is no moment around z axis or x axis by 
symmetry. The moment about a fixed point O 1 on the z axis is 

+C S [15a] 

It is not necessary to determine when the force f~ and the moment 
77?f are the desired information in dynamic problems. For a sphere, 
all the pressure forces acting on the shell passes through its center, 
hence produces no moment about it, i. e. , 

~ *z + C s - 0 Cl 5b] 

Therefore the moment about O’ is simply 

Cl 5c] 

This statement can be easily shown by integration of the moments due 
to pressure on the wall. 

E . Pitching Oscillation 

Consider a pitching oscillation of amplitude ^ around an axis 
which is parallel to y axis and at a vertical distance J? below y axis 
(Fig. 3b). In Figure 3b, it is clear that 
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Z + J? 

COS 7\ — ■ 

JTzbJFTx^ 


sm > = ^ 

J(z+J) 2 ^ x z 


The radial distance of any point (x,y,z) from the axis of ro- 
tation is tJfc+jf+X* • The velocity components on the sphere due to 
rotation are: 

U s = ( rJ{*+J?f+ Qtf) cos A — (z+J>)6y 

W s ==—(-/ (z-t-tf+z* 9j) s/n A = ~XOy 
The boundary condition on the wetted sphere, R, is 



— [ (Z+J?) C*S& — Stnfp c*S& ^ fyl f 


= &yj ( £££ 0) 

This is equivalent to a translational oscillation of amplitude 
in the direction of x. Since the boundary condition on the free surface 
is the same in the presence of a fixed gravitational field, the result for 
translational oscillation can be applied. There is an additional static 
tipping force which can be obtained by integrating the additional static 
pressure f> over R (Ref. 14) 


/ 


p - ps x % 

[l6a] 

F s = Oy f>$ (ft cos(n,x) dS 


= %P9 fff c/, ',(?*) dV = $f>3 % - 
% 

Cl 6b] 


This force acts along an x axis rotating with the tank. 
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Similarly, there is an additional moment 

M s — / p cos (r>,*) ( Z+J>) c/s -J p cos (n, z) X c/s 
R ft 

= f>9^ ra*Jt 

The total force along an axis x rotating with the tank is 

6 - F* 

The moment about O’ is approximately 


Cl 6 c] 


fl6d] 


n?j = F* Jt +m' = F X J? [l6e] 

The total force in the horizontal direction is still . When there is 
tank fixed axial acceleration, the method of superposition presented in 
Reference 14 can be used to determine the x-force. 

An equivalent mechanical model for sloshing in spherical tank 
is given in Reference 15, but unfortunately the extrapolation to include 
damping was not as successful as in the case of a cylindrical tank 
(Ref, 2) and could only be used for order of magnitude estimates (Ref. 15), 
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NUMERICAL METHOD 

A. Approximate Determina tion of the Kernel F unction at a 

F in ite Number of Points 

Numerical quadrature formula will be used to replace the in- 
tegral Equation [6] by a matrix equation. There is a minor difficulty 
due to the presence of logarithm’s singularity at P — I or ^ , 

the latter of which is the integration variable. In the original manu- 
script, an attempt was made to devise a more sophisticated quadra- 
ture formula, expecting higher accuracy. Unfortunately, it seems to 
contain integrals difficult to express in known functions , or require 
very careful process of taking limit under the integral signs. Further, 
the apparent higher order terms may be actually very large and not 
negligible. To reduce total effort, the present numerical scheme 
based on midpoint formula is devised. 

The integrals are divided into N equal parts (N= ZO will be 
used) and the field point is one of the centers of the intervals. A 
simple midpoint formula will not be applicable when the logarithmic 
singularity appears at the midpoint, but if the interval is subdivided 
into four intervals (or more) the error may become acceptable. For 
example, consider the integral 



= A f'fj ~ * = - /• 693 /S <4 


[ 17 ] 
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With four equal subintervals, the midpoint formula yields 

S = [ <2 X (s~ A ) ->-£&> (-§■*)] it = [aJ„ a +( £/,, 3 -/>, #)a] 

. [18] 

= A Jfrt £> — 1.5301 A 

The error is .16^ . For A = 1/20 (N = 20, b=-l), the relative error 
is less than 0.18%. 

Let b ~ 1 , 


h, l ft ' ft) d ft ft 


[19a] 


ft- Sfifs 


fl9b] 


ft~ ii =J fih *' ( ’fi.fi) 


(note the order of i & j) 


Cl 9c] 


f20] 


then Equation T6] can be rewritten as 

lb ft 'rrf.fr fy'is 

Two similar numerical schemes will be presented. The first scheme 
was actually employed in the example, while the second scheme is the 
alternative scheme requiring much less computer time. In the first 
scheme one evaluates and //■* at N x 4N points, assuming four 

point midpoint formula, i^ point on the free surface is located at the 
midpoint of the interval (i= 1,2, ... f N). j represents the inte- 


gration variable located at the midpoint of the subinterval (j=l,2, . . . , 4N), 



Thus 


16 


P: 


£N 


~N 


- 4 +(<'-*) * 


U? = l' [21a] 


p / . (±r! ^ . / * , \ ^ 

Q SN * -?,v ~~ s +( t * 


(L*L ; =/2 [21b] 

To describe the second (alternate) scheme, consider the whole 


square domain to be composed of N x N square subdomains. In all the 

Jo> fo) 

diagonal squares, and M , are evaluated as in the first scheme 

(four values in each square), but they will take the value of the functions 
at the center of each square in off-diagonal domain, which are also eval- 


uated. These data will be denoted 




by ^ ■ h ' 


M 

V 


Total number of 


(o) (o) 

evaluation (both ^ . and M. ) are N x 4N in the first scheme, but 
N x N + 3N in the second scheme (4 point midpoint formula for diagonal 
integral). For instance, consider 

In the first scheme, a, p. being given by Equations fziaj, f2 lbj , 

>' 1 ' '& 


M, 




v r 
= * I 

J *A~ 




(•) 


t & H ‘i 




= j A. ? //'" 


— If A -T ,. (0) — 

^ ^ 4p/ //j.* 

>»/ * ** <t j~( f V 


[ 22 ] 
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In the second scheme, first let = Ja , ft , ft , ft given by Equation 


£2 la] 

N f (») // 


[23 0 ] 

S. . , = 2 ^ X fJ* 

*-' J i-*-*.*-i.**i-.**k * •>* •* 


[2 3a] 

^ A (*) 

J*! ftj 

, i-A 

[23b] 

C * u (b) y _c_ ~t 

™ *~£ Mf- * r 

s 

[23c] 

\<>= 4 ft ^ 

/ <- y jP & 

[23d] 


Equations [23a] — [23d] can be condensed into the single formula with 
ft given by Equation [21b] and one finds 


,t 1! f i >?» 

j-, 4 %■ u y 


f 23 ] 


Similarly, the second integral on the right side of Equation [20] is 
I 0 *# 5 , c <-4 


[24] 


where 




[24a] 
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4t 


. A. 46 (0) 

f ± Z 4 

4 *ji ~ 4 %' 


Jo) 


*(•) 


to) ^ {o) 

In the second scheme, replace Ct. by and H,. by f/.. 

(t f i J 

The integral Equation [6] is therefore reduced approximately to the 
matrix equation, 

f[C]=-fM] -[C] [D] 

the solution of which is 

Cl] is a unit matrix. [C], [Mj, [D] are square matrices of which the 
elements of i^ 1 row and l** 1 column are C.^ , M ; j , D.^ , respectively. 
As a check of accuracy, the symmetric property C-- = Cy should 

9 * 

hold approximately. Then we can use the average value for the cor- 
rective term in the kernel function, i. e. , 


[24bJ 


[25] 


[25a] 


W} 4 ( r'.f/> s <y=?(<v+9J 


The kernel function is therefore 


[26J 




y 

C t y being known at discrete points corresponding to both 

given by equation £zia]. The difficulty of the problem, however, lies 

(0) 

in the accurate and rapid evaluation of the function J J u and ft . . 

* V 


C 27 ] 


19 


at' 

B. Determination of Eigen Vectors, / and Natural Frequencies 

The eigen function takes the value on F, which is governed 

by 

f t> 

where 



a£a_ 

3 


, 3 


being effective gravitational acceleration 


Let /”(f>) %(f>) .then 

f\> n<r-r>l /V »?' 

Analogous to Section A, the matrix approximation of Equation [ 29 J is 



where the factor 1/2 on the right-hand side is in agreement with Refer- 
ence 1, since the strength of the Green’s function has not been doubled 
in this paper. The elements of the matrix A is 


~ tti-k + C (k 


t ¥* & 


A V* / ^ U) \S\ 

H i' +c * =&k +c * 




where i, k corresponds to ft , both given by Equation |^21aJ and 


[28] 

[28a] 

[29] 

[30] 

[30a] 

(30b) 


both vary from 1 to N. 
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In the first scheme, is not evaluated at the center of any square 

subdomain to reduce computing time and is approximated by 

H U + = 1 *> N) f30c] 


where 




(property of symmetry) 


(o) 


In the second scheme, replace //,. by h (i = l to N, j = 1 to 4N). 


The largest eigen values 


v ' V 

N 




of Equation 1^3 0] yields the 


least resonant frequency parameters. Sl n . And the eigen vectors will 
be employed in evaluation of the force response. 


C. 


Evaluation of Force 


The sloshing force for translational oscillation is 




u -fi u ? a ^ - 


where 


[31] 


^ = U**+**'*f ~ 7 f3 ; m* = T- U+3% - 4] 

°(„ = <rrb z J | /”((>) Jf>= 7t~f = vrb 2 .Xj fft 2 ■ A 


A -$„{-&) fa j ftp) j^ 2 Jf> = ?rb l p* = 7r6*J2»(£) ,2, fft ft* 


[31a] 


[31b] 


[31c] 


The nondimensional force 


p* — coo , tfOo a*) 

a " 3 S »•' su(Jg-t) 

D. Precision Problem 

The functions , H g have been first expressed in terms of 
complete elliptical integrals of the first kind, the second kind, and 
the third kind and of simple elementary functions (Appendix I, II). 
The elliptical integrals of the third kind are expressed in terms of 
Heuman's Lamda function A 0 (Ref. 12), which is again expressed 
either in a series form or in a close form of incomplete and com- 
plete elliptical integrals of the first and second kind, i. e. , 

-%r[ e(*)F(\ 8,/0 FK (t)E(fA') -K(t)F(p,&')] 

In ^ and , a serious precision problem occurs due to almost 

complete loss of significant figures in subtractions for ft, ^ both 
small. At first, the series form of the Lamda function was used, 
but it was found that the series is very slowly convergent when the 
parameter is near unity, especially if double precision or twelve 
significant figures are sought. Then it is resorted to the iterative 
methods for evaluating elliptic integrals (Ref. 13), which converges 
to 10 7 within four or five iterations. Although the complete elliptic 
integrals can be computed very rapidly, the subroutine NEFF (Ap- 
pendix III) for incomplete integrals and a difference related to it 
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consumes 8 seconds (used twice), while the total time for evaluating 
7 ^ , H # is only 25 seconds at each point, all on the GE 225 computer. 
Longer time would be required for higher precision as the number of 
iterations increases. 

To increase the precision, analytic subtractions are made so 
that no significant subtraction remains, if possible. Noniterative sub- 
tractions in which four or less figures are lost are acceptable if four 
or more significant figures out of eight (single regular precision on 
the machine) is desired. The technique can be illustrated by the 
following cases: 

(1) Let (A — B), the difference of A and B is small but can be 
expressed analytically without subtraction. Then, for ex- 

e -g-> A=2 , 6 = 2+Z , S «/ , (A-B)=Z 

(2) Let K n ' s be small (positive) quantities containing no sub- 

traction, then — / should be evaluated 

by repeated application of the simple relation that 

(/ + fi )( f +£ 2 )-/ — ]£ f +&2 f 

(3) To subtract a desired quantity from a known function may 
require a new subroutine for this function performing sig- 
nificant subtraction analytically, e, g., NEFF. (Appendix III) 


should be evaluated from 


zL (dz*l 
Ja/b JX+J& 


ample '^ - ^ 
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Aside from relatively mechanical operations, the device of DKEF and 
NEFF subroutines, the following relation was expedient (Appendix IV), 

(±± SUji jk(*> I 

^ ~ * 1 1 ( /-«') f r*VI 

JT n and Ji /( j are defined by Equations £lV- 6 j, pV- 8 ], respectively. 

It is noted that, after a small manipulation, direct numerical 
integration of the integrals F$ , // y at sampling points of the entire 
domain of p £/ ^ was also computed by Weddle's rule. Although four 
or more significant figures can be obtained, it is deemed too slow over 
the major part of the domain. For instance it took about 5 and 2-1/2 
minutes respectively for F 3 and H tj on a GE 225 computer with 384 
intervals, or a relative error of about 10 “^ otherwise at a point near 
the right lower corner of the domain ( p. , near unity). These values 
at sampling points are valuable as they serve as a good check on the 
present computer program, which evaluates F j , F £, F 3 , Hqi» Hq 2 » 

Hqj, all together at a rate of 25 seconds per net point (on the same 
computer with an accuracy of four or more significant figures). 


1 
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EXAMPLE: FUEL SLOSHING IN A QUARTER -FULL TANK 


First, and H^are generated, then the matrix Equation is 

solved. The corrective part Cy to the kernel function obtained is sym- 
metric almost to four figures (Table I). The relative errors in the sample 
points are less than 0.3% or better. Since these values are quite repre- 
sentative, the values of Q. at other points are not shown in the table. 

Next, the eigen values and eigen vectors of Equation [30] and then 
the force response of Equation [3 1] are calculated. The calculated first 
four eigen values are 9.48863, 2.0591201, 1.2003387, 0.84773955, re- 
spectively. The corresponding frequency parameters are compared with 
experiments in Figures 4a and 4b. It seems that the values are well within 
possible experimental error, although it may be slightly less than the 
actual value, noting that natural frequencies are somewhat smaller for 
larger amplitudes of oscillation. 

The constants needed to calculate the force response are compared 
with graphical values given by Reference 1 in Table II. Since the coefficient 
is in agreement with Budiansky’s value, the main difference lies 
in the value of first natural frequency for frequency range in its neighbor- 
hood. Since graphically interpolated value is less reliable, which is also 
confirmed experimentally in this case, only the present theory is compared 
with experiments (Ref. 15) in Figure 5. The difference between theory and 
experiments, perhaps, is essentially due to finite amplitude effect. But 
the agreement seems to be quite reasonable. 






TABLE II: Comparison of Constants with Data 

from Reference I 






CONCLUSIONS AND DISCUSSIONS 


The i" 6 s 6 Tit theory 3.11(5. coiuputGr program seem to yie 


f y l p 1 rl 


satisfactory predictions of natural frequencies and force response 
in comparison with experiments for a quarter-full spherical tank. 
The computer program is expected to be applicable to other liquid 
depths, although not beyond improvement in efficiency. The results 
also confirm the theory that the kernel function is related to the 
Neumann function on the boundary and that this function can be con- 
structed by adding a corrective part to a known Green’s function 
numerically for practical applications. Extensions to other prob- 
lems may be possible, but one must resolve the precision problem 
if it exists and one may also find a more sophisticated numerical 
scheme to be more desirable, either in accuracy or in efficiency. 
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a = 



b - 
d = 


DC9) *• 
DKEF = 

F = 
^ = 

^ - 

z(tr> = 

<? = 

A. - 
“V 

Grea) = 

G//5Q) = 
4(7? <3) - 


- 

F(/?fO = 
F.<FF> = 


NOMENCLATURE 

radius of the spherical tank 

area of undisturbed free surface 

area of wetted surface of sphere 

maximum value of p , radius of free surface 

tank diameter, 2a 

(e(Q)-K(9)) 

q 2 (c.f. Appendix III) 

a function in the computer program (c.f. Appendix III) 
the undisturbed free surface 

horizontal force acting on the tank due to fuel sloshing 

___ ( 0 ) 

- % ( note the order of in ify ) 

integrated kernel function related to 
effective gravitational acceleration 

dfrfc h (ft, fy) 

Green's function of the second kind for the spherical bowl 

Green's function of the second kind for a sphere 

additional part of Green’s function for spherical tank other 
than half -full 

integrated kernel function related to ft(P'Q) 
integrated kernel function related to G(P,Q) 
integrated kernel function related to Q+(P, Q) 

•fWfi H(fi.fj) 
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K(S\E(S) -r 

K = 

n = 


NEFF — 


m ) = 

p(i<p) — 


point of integration, except I being the unit matrix 

complete elliptic integrals of first and second kind, 
respectively 

total mass of liquid (fuel) 
outer normal 

a function in the computer program (c . f. Appendix III) 
( Km- f ) (c.f. Appendix III) 

a ring corresponding to p( O <P, &) 


% 

% 

a,i 

Qs 

r^e 

R 


- iJRL 

p+p' 


V 

V (pp - Rf-t- a* 

analogous to P but related to Q and I, respectively- 


defined by Equation fl-lOa] 


spherical coordinates 


=r the wetted spherical surface before sloshing unless defined by [l-2 


R(RP) 

R Ru R.. 
ta J ** > 7 

dS 

dS 

U 

X 

X F 

°/7 


defined by Equations [l-10hj, [i-lObJ, (ll-5b}, respectively 
distance between the points P and Q 
element of surface 

dS/ dQ -+ f>df> on F 

horizontal displacement of container in the x-direction 

r cos© 

vertical distance of free surface from center of sphere; 
positive upward 


defined by Equation [l2b] , £ <£(l)c/S x 


Pn 

cos y 

A, , A, 

Sp, i) 

Tf« S) 

p 



P 

% 

6 


t 

UJ 

*>n 


defined by Equation [l2a] 

-r angle between the vectors OP and Ofy 

= [- ± J zg , respectively 

ss Heuman’s lambda function (Ref. 9) 

— C Oy 

3 

= complete elliptic integrals of the third kind (Ref. 9) 

= radial distance from a point on the free surface to the 
center of the free surface 

=■ f> of integration variable 

density of liquid (fuel) 

= velocity potential, 7 = j? J 9^ being the velocity vector 
= nth eigen function 

=* nth integrated eigen function related to ^6 

= jf %<t> 

— frequency of oscillation 
= nth resonant frequency 

— tOrfa/q } n th resonant frequency parameter 


Subscripts 


F related to surface F 

i , j , k related to ft , , ft , respectively 

I related to integration variables 

P related to the point P(r, 0) or P(f > p . ) 

Q related to the point or Q(fa) 


R 


related to surface R 
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APPENDIX I. ANALYTIC EXPRESSION FOR N a ( 


The Green's function of the second kind for a whole sphere 
(Ref. 16) is 


where 


R — <J r*+r'*-2rr'cas? , A? = Jr' V -^1 - £T 

cos/ = cos </> cos<p -hsio <P <j> eas(& - e') 


[1-1] 

[ 1 - 2 ] 

[1-3] 


When P and P ! both on F, 


GAR,P)- ' 


-jff 7 ^ 

J{>*+P' z -2f>(>'c*s(9-9') J f> > 2 +f)Zf -2 pfa cos {&-&')+£ 


+ <* 2 2 

S-l £ -r* A /» 




-/,')/ A 4 " I 


[1-4] 


* " z f -ff* c os(e-e') + J -2f>pt^cos (e-e'HB 
Using Equation [l-4j , with q- — q~ $' 

_ _ 

H 0 (Pj P') — f G a (p,F>) eosQ'Jo = / G.(r,r'; cos<TjtJ> <p') c*s<r- e/g- 


- H ot + H 9Z +H t 


03 


[1-5] 


Making use of a new variable yS Jp 


H c 


27T 


cos q- do~ 


° rj f>\ p /J -a pp'cos (&-& ) 


and a new parameter 





fl-6] 
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where 


2* 

-2 


m) -f , = = 

£(*) = 9* srfp <*f 

This can be evaluated more accurately by 


/%,= 


_ / 


r S, [~K(% )-zo(%)] 


where D(9,) is DKEF (3) given by Appendix 1 1 1. 
By taking limiting process. 


= H ot (f>,o) =o 


Similarly 


27T 


Hoz ~ 47r 


a cos o~ dcr 


0 J ( Zpp'cf-b 2 f>f>a*cos <r 




where 


?,= 


:a dpp' 


$ / 


J (pp - ~i*f+ 

Also, this can be evaluated more accurately by 


H °* srtjppr I ZD(9 ‘}\ 


[l-6a] 


[l-6b] 


[l-7] 


[l-7a] 
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When p or p is zero 

HCo.f) =H 0i (f' °> = o 


From integration by parts, 


[l-7b] 


C'P pp w<r + J( pf>'—b 2 ?+2pf>Q 3 (/-c,os<r-)-t-Zp(f>-p^ cos<r ^ <r 

/ _ d<s" 

^ na L V (fp'-lff-hzppcfQ- ~coso) +*p(f>-p'f P'P')*] 

7f 

= Tza~f o / [fpJtfp-lfftzppa^O- w<r) -hig(p-f? + ppt 2 2 ] • 


I 


C J~ ( Pp'-b*T+zpf'<£( T- COS <r)+z*(p-f'f -f+pptoso] -stoV* 

♦ 

(-p 2 p f2 )(ctsT-\ )(cas f?p*bff+Jj?({+p'f— 2 f>p'a*c*s <r 

pp' /’ 7r I (ftp ') 2 — 2 j> p d cos oj— G?i> + (ffp casC~j sJr?Q 
W 0 >) p*p' 2 +b 4 -+‘Z^( P+p'f-2-fp'ci COS cr (~P* f> /Z ) (cos<r-),)(casa--)i z ) 


do- 


dr + 


- ,* * /* , 


+ -pC( — ±££Jsss Jt 

zfra Jo -py*( cos <r-x) (cos<r--A 2 ) 


= I.+I* +1* +1+ 


[l-8] 
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where 


Sir? Z <J~ 


-/+ 




( cos<j '-) l ) (cos <r-)i 2 ) (cos it-- Ji,)(cosO'-?ij) ( cosu~ - A,) (cas(T-^ ) 


27r« ^ (-py,*) y ( p+fff—zppc? cos<r~ 


■ (-pp.) j fi- 


IL 

,£ 


tfo ~+h*+ z£(p*+f>' 2 )-2C?l 2 ] 


«*( rYJ. ‘‘ 


c/t' 



* rtrlfr { iTjpp) w+e+*<r*f>-«re J *rv f 

[l-8aj 


By using partial fractions, the variable jB and the definition of #"(**$?,) 
(c.f„ Ref 0 9 ), one finds . 


j 
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fr / 

T f cffl-rfiacK <r r / ; -(A^)<«so-f^,A» ] 

fy 3 +lf+^(p+pf-zf/Sa COBS' [ (cos<r-\ ) Ccos^-A,) (c*s<r-\) (casr-X^y*^ 

_,_i tfAb* +^(p+f j- <?&*!- f , Pa t a>s<x- ( i / v 

* m fTj 0 4 f>*p*+b 4 -t z?(ptf>'f-2f>{**ax<r ( VA* *““»*-*/ w-*«' 


Aj^)* / Af ^2 \ ) 

A, -A* ( *«<r- A, w-A* /} * 





j ^ X ' _ Aa \ 

^p'\b 4 -h^f(p+P'f-Zf>p'o^cos(A C< x ,r -*I ox<r-\J 


f+ 


- -¥-- r o.-w+.-rV - — \ it *■ 


X, 


A ; — A z 


c*scr- A/ 


--^rr+s+*w-«vi[^ jfftjfe . i)+ 


U ^f-v/ ' } *)] + ifrf^] 


+ ^ «*> *0} 


where the well-known complete elliptic integral of the third kind is 
^defined by 


zf 1 , at*; 

j](^2)= f ^ — = / — 

7o (>-<** stn z fi)J / - gsm*p Jo ( 


'* sr>U ) 


Cl -8b] 


[l-9] 
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z * ~ F y*)l-(o-6*+f>pc«<r) c/<r 

~ Tfta. (~FF^ ^ ft] 


[l- 8c ] 


j / / 7 r / r (q-b^+p/oosT + 0 a 1 -£ L + PP c os <r)(\ ~ (\+\t) cost ) j 

* 27 ra ff > ' Jo ft'* l(CK 0 ~-X,)(c 4 S<r-Xj) (cAs<r-X,) ( cm cr-^i) J <r 

_ f ,7t 

2TT4 pt>' X l~ (X ' +X * Vf + (c*sir-x l )(<iisa‘-X M ) [ff (*+ >] + 

■+ (X I ^X Z ) X t Xzpp (Q^-l?) J ^ {/<f 


zJq I ( VA») 7t ~~pft Jfti r /»A | 


[l-8dj 


where 


Hence 


*■ 








[1-8'] 


+ i»p- "■ y A+W' r - tpp'w^-jtp Jfci Cff *,+■£] | 


I 
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When ^ or ^ is zero, 

It was found that there is a precision problem in Equation(j-8i| 

✓ 

when f and ^ are near zero. This might be anticipated as there is 
a very small denominator proportional to 07° and the result is 

expected to be small in view of Equation fl-8'j. After somewhat labor- 
ious manipulations with Equation to resolve the precision problem, 

is obtained in the following form (with 1): 




*** f pp'Qs - 

* 2 ^ ^ - } -y ’i 


[1-103 


whe re 


<5, 




ZCtJff 7 


[i-lOa] 




(l-10b] 


* 


M 

I A 2 


r £f<fc) +&&)]+ 


K(9 X ) Ol-O 

** a,-**n z * 


[l-10cj 
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R h . = R, +R Z + +r ?4 


fl-lOdJ 




-*» W^y<jrV>J -f 0*-*r>? } 


[I- lOe] 


/?, = .-3L- 

2 a, 2 


tt / *(/>%r*)+rr * p%e'\ 

** = *pp' | Ft*, +*f * j ' 


Ci-iof] 


[l-10g] 


*V«= y? 7 ? 7 ?v^j 


( I- ! Oh] 




<7T 


V^rOa*) 


(l-10i] 


A, > / , >i a < 0 


fl-10j] 


It must be noted that C\~0 , (£"/), (/f^-/ ) are evaluated not 

by direct subtractions, but by accurate formulas. 


( R il -0 =f>Y 2 + aft 4-$(f+ff 


[i-lOkJ 
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[I- 101] 


I r-r'\ 

rr' + rp' + **a 


fl-lOrn] 


(£-/) = 



JIzJ+2R^ Hf+f'*) J / (f>Y z ) 



fl-lOn] 


(p-f') can be calculated without loss of significant figure. It is 
also noted that for small Xp > one should replace (t^-h z a") = k*(fi-<£) 
by —Zpb or ( t—G 2 ) by - with b = f w 

It is recalled that Budiansky's technique of differentiation 
under integral sign does not seem to lead to simple results, due to 
the presence of non-zero , the relative depth measured from the 

center of the spherical tank. 
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APPENDIX II. ANALYTIC EXPRESSION FOR -f( 

The outer normal derivative of Gf,(PP) on the free surface is 

# - - ifr I ~^ (z ’ *■ p (t *- f> z > + 7^^) [>"'■=»*-'-*< J + 

U^'+rr'**)- <VVA«V*Wjf 

When both P and P’ on F , 




-f — T ( a*r + rr' 2 ) - ( Sr'-t-A*/"') cos'}] \ 

ax'r*r' i (c*s i 7-,) L * 


where 


r*=p z +z F 2 ,r'*=f>'*+4 

For P, P ! both on F, 


3(pyp')= ( 4^1 </<r 

7 0 dx » x=st r y =r^ 

= 3 + 3* + 3 3 

Using the same technique as in Appendix I, one finds 

/ 27t 3 

3 (f> ' ^ "J o -jfc' ^ ( 7"' f?) c * s<r J<T ~ 

- 7 7 ^ I ( f; 


In- 1] 


[II -2] 


[n-3] 


[II-4] 


where <2 3 — r*« = ct(tf-p 2 ) 
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It is noted that a special case of the reduction formula in Reference 13 
can be used to evaluate the following integral, which occurred in 

it _ 

r I-Qtsin'ff* £(v 

Or, use 

4* 

X c-K x ^*f A 1 J * 

which can be checked easily by differentiation. 

When f or p is zero, 

3,(0^)= 5,(f,o) = 0 

To increase accuracy in numerical evaluation, Equation [ll-4j is 
replaced by (with b = 1 ) 



+*<**•>] <t+f»o-r> 

/ / ^ ^ 

0~f > p-p and b-f>f> can be evaluated accurately for known 

discrete values of p, p and is evaluated by 


J 4f>p'a* + U> 2 -pf>'f+ $ (f> -/>* 


Next, will be expressed in terms of elementary functions. 


[ll-4a] 


[II -4b] 


[ll-4c] 


| 

I 
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iW 


Z(f'f> = 


-/ 


47ta Jo r*r'* ( cos z Y-t) 


~ [ casft— r"i cAS CT c/tf- 


7T 


-=Sl{ 

J O ^ Si* \ /'*« 2 J 


( ^ 4 f>f>'oos<T) Z -( *?+{>*) (z^^ff 2 ) 


_ ~ Z F f 7T 

£7ta 


ti± z£l 2*£± M 


[— + r tii zfH cosa ~+ &£ c i 

l Pf Jo f*r* [ cos <r--\j [ 


cos <f'-) 3 J 


’-Jr\ 


_ — Zp _ r c 

■ z ‘/ > / 0 ^ Za f , y 


t r jJg£^£V % , /.-Vrtl— i- - 

v<vv r <-/?'> f -W 


— * r 

2 pya <*,->,) L { -rr ) F J V^ : 7 


-ZlB f _J_ r t&Jzjj +*£(&*)+&* n 

z ff > '°- < 2 l zg+sgCf+^+pY* J 


+ 


—L L £ 5 e*£) J4+(p z +p'*)4 +&>* . 1 1 

i$+*?<r+r*> + fr x JJ 


Some further manipulation is required to avoid precision problems 
for ^ small or ^ small. One finds 

7 -*• Uf E£-fJJ± L £ t*k ) CfrV*H ^Vll] 

* / Z L \ Z f\ J 

. j - . z 2±£_ . [ii£l ) 

zp+f+fk \(rr > I i^i ^ i\* F \f?tj +z//-rf fa ]j 

= 'JTjf+Zi %a Ty+p- 


[II -5] 

[il- 5aj 
fll -5b] 
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When p or p — > o 
? z ( 0 / f)= 7 jf>. 0)=:0 

Finally, will be expressed in closed form as follows: 


[iI-Sc] 


= -. 3 - F [ -( r^a 1 ) ( V +f>p'cos <r) ] . 

ZK ° J° J -hzpf'^+z^p-p'f-zppa 2 oxer 


cos o~ </<r 


r -Zf(p a +fr z - 2pf ' cj>*<r) - f 3 f' 2 (/-ca$‘cr) ] 


.J£ 


7ra L z *ffF a I J/-? 3 ^p} j[ ( p e - ] + p*r'{x,-x)(u*<r-x)’ 

■ [\ [frVa ^ 1 + r*(a+r' 2 ))- (pf(f> z +p ,!t )+p*p'*) ] + 


ZlL 9 * 

Z 7 ra 


+ f>Y* ( *r^(c~<*<r -*) fa ( (r *+ °?>*F + r*(e 2 +r' z ))- tjjrfe 

7 = = • 

a <ff>y « ^ jz^+^^+p' 1 ) +p* r * 

• ~ ‘ [\/( Jjc+p 2 +a)*F -t-(zf+p l )(a*+jg +P' 2 )) - 

-ry,^, *»%?* ]. TTj^y — ± = ■ 


Ufa , a) j 
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_ Z F % 


zita 


5^5* pWM <(%)+{[/- 


(_ pp - ^-y ^ +^(p 3 +r^ +fY 2 

z? +Zr*(p 2 +p' 3 )+fZp'* 


¥('f+p? [(-^+J4+&f**f~>+f ll p t X* J £+$ ^VrV^V^ - 
-Cy^/>V(*V^>VV^- Tf(77x,^') - 


_ j. r, , f PP- 4 +*?(?*+ r^PY 2 -, / . 

2 1 ’Z+^Cf'+r')*??' J ^ w>* 

•f(V +l$+*£(f*+f*)+f>*r*‘ )(■ 2 *r iz F (*«*+*f*+{*)+f(e*+r*)) + 


+ (*F+f z +* 2 ) ( *£<p*+p*) +/>y*)] If(-j^j s f 11 " 6 1 

There is a serious precision problem for ^ ^ small in 
Equation £lI-6]. After manipulations, the precision problem is re- 
solved by employing the following equivalent form. 


7 Z f q - 

s 2Tra l (f>fy v * 



Cii-tI 


where 


F„ = -\ r*W«‘;7£-f £ [B,+2(*?*p’+a‘)B,*B l B s ](K(),)+Tl)\ 


fn-7aj 


48 


A- 


(f s -f' z ) 


( ±tn 

Fit « CXta+ZF+fV 


5* = 



z p F/%a 



-Jt 

z hr*+«* 



Where the terms in the inner bracket could be replaced by 


*/ +#** 
z pFp 2 -+f?t* 


for higher precision, which seems unnecessary as the error in B z is 


2 * 

sufficiently small in the critical range due to the factor p' - f>* 






o s { -^r -<r~r‘> 




ri£ln£i 


< P+P') (p*+ Zp) R ia Z P*+Zp +f*ta] 


1 J 


(r-r) (4 t£LL : £i£±rll +(*?+*?) ail 2 

P'( *?+ p 2 ) Kn (*F+P 2 ) +PY 

<Xz£l ) 

P«ta * 


[ll-7b] 


Dl-7c] 


[Il-7d] 


[Il-7e] 
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where further manipulation may be needed for very small * F to 
avoid precision problems in the domain of small f and (*' , 



-(f-n ( z ?+f 2 > [- %*+ 


(z*+p 3 +a')(z?£-hf) _ 


When jo or f>'-+ 0, ^ — ► 0. 

It is important to note that whether f> or the sum of 

and always approaches zero as 0. Therefore, 

for a half -full tank H=H 0 which is in agreement with Budiansky's 
kernel function aside from an apparent factor of two difference 


r«-7f] 


fn-7g] 


mentioned previously. 
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APPENDIX III. SUBROUTINES DKEF AND NEFF (WIZ PROGRAM) 

DKEF = DKEF (k,k, 1,1,1) 

NEFF = NEFF ( ^ , TC , 1 , 1 ) 

The unity arguments are actually dummies, while the five argu- 
ments represent five outputs. For DKEF, the outputs are K(k) » DKEF (1), 
E(k)= DKEF (2), (E(k)-K(k) )/k 2 - DKEF (3), K(k)-*r/2 * DKEF (4), and the 
number of iterations = DKEF (5). DKEF (3) is not obtained from DKEF (1) — 
DKEF (2) but is obtained after a significant analytic subtraction in the pro- 
gram. For NEFF the outputs are F( k ~ fc ) — NEFF (1), E ( y$ , k = tr ) 

= NEFF (Z), (E( p , k) - k sin p ) =- NEFF (3), the number of iterations for 
evaluating F ( ft ,k)- NEFF (4), the number of iterations for evaluating 
NEFF (SJ^NEFF (5). NEFF (3) is evaluated after a significant analytic 
subtraction in the program while NEFF (Z) is simply obtained from NEFF(3) 
k sinf* . Although k' — J 1 -k^ does not appear in the functions sought, 
it is calculated from a formula without subtraction, as one can easily see 
significant figures of k’ would be lost if k is near unity. The basic formulae 
are all given in Reference 13. For complete elliptic integrals, the iterative 
method based on geometric and arithmetic means was employed. For in- 
complete elliptical integrals, the iterative method based on inverse order 
of transformation was employed in order to construct NEFF (3). The pro- 
grams are written in "WIZ 11 language for GE ZZ5 computers, which is 
analogous to "FORTRAN" for IBM computers, and are given on the following 


pages : 


WIZ SOURCE PROGRAM 


SEQ LABL TY 


STATEMENT 


C ZE NZE PL Ml 


4oo_dkef_. 

U 01 _ 

402_ 

403_ 

4o4_ 

4lO_* 

4 1 1 _ 

4l2_ 

420 _ 

42 1_ 

422_ 

43°_ 

43 1 _ 

44o_ 

44l_ 

442_ 

443_ 

444_ 

445 _ 


_ARGjft>KEr{1) $DKEF(K,KP,1 ,1,1) _ 

_ARGPi»KEF(2) 

_VA#1 ,VB#ARGP,PI # 3 . 1415926536 

_QROO#ARG*ARG/( ( 1 &ARGP)* ( 1 JcARGP ) ) ,KN#QROD .N0#0_ 
_KK#0. 5*vROO , SUN#-0.5* ( 1 &KK) 

_KPjfVB/VA,SUM#SUN,PN# 2*SQRT. ( VB*VA)/(VA&VB) _ 
_KN#KN*KN/( ( 1&PN)*( 1&PN) ) ,ROD2#KN ,ROD1 #JROO 
_QROO#ROO! &ROD2&ROD1 *R0D2 
_KK#KK*0.5*KN ,SUN#SUM-0.5*KK 
_VAT^0.5* (VA&VB) »VBT#SQRT. ( VA*VB) 

J/A/fVAT , VB #VBT , NO #NO&1 
_ABS. ( (QROO-ROOI )/QROO) -DELTA 
_ABS. ( (SUN-SUM)/SUN)-DELTA 
_KNPH#QROO*PI*0 . 5 ,FK#KNPH&0.5*PI 
_DKEF(1)#FK $K(K) 

_OKEF ( 3 ) #SUN*FK $( (E(K)-K(K) )/K/K )_ 

_DKEF ( 2 ) jj*FK&ARG # ARG*SUN*FK $E(K) 

_OKEF(4)#KNPH $K(K)-Pl/2 

_OKEF(5)#NO 


ANY 


52 


SEQ LABI rr 


WIZ SOURCE PROGRAM 
STATEMENT 


C ZE NZE PL Ml ANY 


500_NEfF_._ 

505_ _ _ 

5'0_ _ _ 

5 1 5_* _ _ 

520 _ 

525_ _ _ 

530J4 

535.J5 __ 
540_ 

5*»5_ _ _ 

550_ _ _ 

555_ _ _ 

55<L 


BETA#NEFF{1) $KAPA NONZERO 

_KAPA#NEFF ( 2 ) , KAPAP^NEFF ( 3 ) $NEFF(4,5) OUMMY 

VA#1 ,VB#KAPAP .BN^BETA, I F #0, NOT #0, TEMP #1 
KPjfVBAA ,PROO#TEMP , I FT#» F 
OCdtOS. ( 2*BN) ,NS#SIN.(2*BN) 

6P#2*BN-ATAN. ( (1 -KP)*NS/( ( 1 -KP)*DC&S SeKP) ) 
BP#2*BN&ATAN.((1-KP) # (-NS)/((1-KP)*DCJc1 SeKP ) ) 
TEMP#PROO/{ 1 SeKP )* ( BP/BN ) 

IFij&ETA'TEMP 

VATjf0.5*(VA4VB) ,VBT#SQRT. (VA*VB) ,N01 vWMOl&l 
VAWAT ,VB#V8T ,8N)feP 
ABS. ( ( IF- I FT)/l F ) -DELTA 
$ IF COMPUTED 


56 o_ 

565_ 

566 _ 

570_ 

575_** _ _ 
58 o_ 

58»_ _ _ 

582 _ 

583_ _ _ 

584_ 

585_ 

586 _ 


DELO^APAP*KAPAP/( 1&KAPA) 
_KX#1 ,MM#1 ,ENKS#(1-KAPA)»IF,FKK#IF # KAPA # KAPA 
_NO2i#0 

_KN#kAPA,DELN#DELO,SSO#SIN. (BETA) ,SSN#SS0 
_SSQ#SSN*SSN , SSP#KN*KN*SSQ , RTKN#SQRT . ( KN ) 

_RT1 yfeORT. ( 1 -SSQ) ,RT2#SQRT. ( 1 -3SP ) ,ENKT#ENKS 
_RT3#SQRT.(0.5*(1&RT1 )) ,RT4#SQRT. (0.5*(1-RT1 ) ) 
_RT5 #SQRT. ( 0. 5 * ( 1 &RT2 ) ) , RT 6 #SQRT . ( 0. 5* ( 1 -RT2 ) ) 
_AC#. 5 * ( 1 &KN ) *OELN*SSQ/( ( RT 3 ART 5 ) * ( RT 1 &RT 2 ) ) 
_AS# - ( 1 &KN) *DELN*SSQ/( (RT4&RT6 ) * ( RT1 4RT2 ) )/2 
_SSH,(«T4.CSH#RT3 

_DSN#AC*SSH&AS*CSH&OELN*SSN/( 1 &RTKN) 


J5 

J5 
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SEQ LABI TY 


WIZ SOURCE PROGRAM 
STATEMENT 


C ZE NZE PL Ml 


5&7_ 

590_ 

59i_ 

592_ 

593_ 

59>+_ 

595_ 

596_ 

597_ 

598_ 

599_ 


KK,fKK*KN ,MM#2*MM 

SSN#SSN&AC*SSH&AS*CSH, N02#N02&I $$ 

OELNijt>ELN*DELN/{ (1 &KN)*( t &KN4c2*RTKN) ) 

KN#2*RTKN/ ( 1 4cKN ) 

vn#-mm*deln/( kk # kn ) 

UN#-MM»DSN/5QRT . { KK ) 

ENK3#ENKT«eFKK*VN-KAPA*UN $ E-K*SI N. (BETA) 

A 8 S. ( ( ENK5-ENKT ) /ENKS ) -DELTA ** 

_NEFF(1)#IF $ F(BETA,K#KAPA) 

NEFF ( 2 ) #ENK3&KAPA*SS0 $ E(BETA,K) 

.NEFF ( 3 ) (jtNK3 , NEFF (4 ) #N01 , NE FF ( 5 ) #N02 
$ENKS IS E(B£TA,K)-K*SIN(BETA),K#kAPA 


10 ENO 


*• 


ANY 


999. 
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APPENDIX IV. DERIVATION OF TT 


N 


For j the formulas 410.01 and 411.01 of Reference 12 


are applicable, in which ~~r — < 0 s ^ & 


1 1 . e. 


7 —7— - 




where 


[IV- 1] 




[IV- la] 


and 


77 .= 

where 


k(*) + $_ ^LLASMlzJj. 

'~* z J Td^F) 


[IV- 2 ] 


. . -/ / 
{5 — sm 


Stt 2 [IV -2a] 

To exploit the possibility of gain in significant figures, JJ Z ^77 * — 

/ ~-r 

K(£) 

will not be computed by simple subtraction of from either 

of the above equations. This difference is defined as 7T n - 

From simple algebraic manipulation of Equations (IV-1) and (IV-2), 
one finds 


K- 


K(t) 


7r w* 


+ 




j (/-«?[/ WV-/] (t-i*')+ 


[IV -3] 


! 


E or p£ , pj small, is near _2L , k is near zero, thus 

Ao(p,£) is near but less than unity. Also for p t - , ^ near unity, k 
is near unity. Equation (IV-3) may still lose too many significant 
figures through A„ -/ or /- . One can further apply the addition 

formula (#153.01, Ref. 12) restricted to the condition that ^tan^» 
tan = 1, i. e. , 

A 0 CPJ) +A. (p. t)= I -h — * '* J * fi 

71 *! OOS*p + tfsin 2 ? 

where ik' 2 = /- jf 2 

Eliminating A.(f,A)-l from (IV-4), (IV-la), (IV-2a), one finds 



K(*\ ; i ( IT «*(-/) (/- 
/ - -| a I- f-* l ,vV<7-oA>(cA-**,> 


AcCk£)+ 


«K(A) (i-iA 2 ) 

(-«*+**) 


The refore 


! 

! 



-or* 




Applying the addition formula (IV-4) again, one finds 


,* / 7r (y ±2. -jKte) | 

J~ «‘( I- **)(**- of) (/-ot 2 ) j 


[IV-4] 


[IV- 5] 


[IV-6J 


PV-7J 


I 
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J~id ~ = (>~ sin P —(A 0 - Sin ft)) 


= ~ ^ , -■ -/w 

/jl-d 2 {!+J I-**) 

[IV- 8] 

1 

o 

H 

>> 

© 


~4~f[ £(*)~F(*)] F(pA')+K(fi) E(p,*')-^-sbp] 


- #■ {** P(*) F(p. *') + PWE(M') + f 0(p, *') * ^_L J 

IV- 9j 

where 


°<*= -*/(-/- A 2 ) , t = 2 e 

[lV-9aj 

D(t) = OQ2 = ( E(*)-F(*))/j? ~ DKEF(3) 

[lV-9b] 

p(t) = PQ2 = ( K(t)~2f) = DKEF(+) 

[lV-9cJ 

D( p,F) = 102 =(E(p,tf)- t' sin ft) = A/EFF(3) 

[IV- 9d] 


There is apparently a gain of significant figures of TT N when 
are small small) if equation (IV-7) is used, provided that the 

first term in the bracket can be evaluated as accurately as the second 
term. This is achieved by employing the subroutines DKEF and NEFF 


for equation (IV-9). 


APPENDIX V. X-FORCE ACTING ON THE TANK BY 
INTEGRATION OF PRESSURE 


Assume a velocity potential 


4 >. — U x + z a„<k,(r,<f>,e) 

v y?=f 

where the first term is a particular solution satisfying the normal 
derivative condition on the sphere. are the eigen functions which 

have no contribution to the normal velocity on the sphere. In order 
to satisfy the free surface condition for sinusoidal oscillations 


[v-l] 




3 ” 


tif * * CO 2 00 . oo , uji 

■ ° r f u * + T £, M T *" ° nF 


one has 


Qrj — 






(*-?)«■ e-'A- 

since i,'s are orthogonal on F, <*n — f dS and /3„- — [x<kds 

J p r 9 J F 

The pressure 


fV-2] 




[V-3] 


The x-force can be obtained by direct integration of pressure 
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since 

-/ t *(**>-*#(]*! 

J Y. \ 



Equations £V-2, -43 are the same results as that obtained in 
Reference 1 through Lagrangians’ equations. 







FIGURE 3b. Moment of Tank in Pitching (Rocking) Oscillation 











DIMENSIONLESS FORCE AMPLITUDE 


z-r /a = -//z 


u/d - 0.0032.0 


FIGURE 


specific 

gravity 

A Water T * / /u = /.0 centipoise 



NATURAL FREQUENCY PARAMETER BASED ON DIAMETER 

9 


Comparison of Force Response for Quarter-Full Tank 
with Experiments by Abramson, et al. , (Ref. 15) 


p ro 



NATURAL FREQUENCY PARAMETER 



FIGURE 6. Natural Frequencies Given by Riley- Trembath (Ref. 17) 




